Theta

# ARIMA-Black-Scholes: Semi-Parametric Risk-Neutral Pricing
# T. Moudiki
# 2025-12-04

library(esgtoolkit)
## Loading required package: ggplot2
## Loading required package: gridExtra
## Loading required package: reshape2
## Loading required package: VineCopula
## Loading required package: randtoolbox
## Loading required package: rngWELL
## This is randtoolbox. For an overview, type 'help("randtoolbox")'.
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: data.table
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:zoo':
## 
##     yearmon, yearqtr
## The following objects are masked from 'package:reshape2':
## 
##     dcast, melt
## 
##  
##  This is version 1.10.0 of esgtoolkit. Starting with 1.0.0, package renamed as: 'esgtoolkit' (lowercase) 
##  
## 
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(ahead)

# =============================================================================
# 1 - SVJD SIMULATION (Physical Measure)
# =============================================================================

set.seed(123)
n <- 100L
h <- 5
freq <- "daily"
r <- 0.05
maturity <- 5
S0 <- 100
mu <- 0.08
sigma <- 0.04

# Simulate under physical measure with stochastic volatility and jumps
sim_GBM <- esgtoolkit::simdiff(n=n, horizon = h, frequency = freq, x0=S0, theta1 = mu, theta2 = sigma)
sim_SVJD <- esgtoolkit::rsvjd(n=n, r0=mu)
sim_Heston <- esgtoolkit::rsvjd(n=n, r0=mu, 
                                lambda = 0,
                                mu_J = 0,
                                sigma_J = 0)

cat("Simulation dimensions:\n")
## Simulation dimensions:
cat("Start:", start(sim_SVJD), "\n")
## Start: 0 1
cat("End:", end(sim_SVJD), "\n")
## End: 5 1
cat("Paths:", ncol(sim_SVJD), "\n")
## Paths: 100
cat("Time steps:", nrow(sim_SVJD), "\n\n")
## Time steps: 1261
par(mfrow=c(1, 3))
# Plot historical (physical measure) paths
esgtoolkit::esgplotbands(sim_GBM, main="GBM Paths under the Physical Measure", 
                         xlab="Time",
                         ylab="Stock prices")
esgtoolkit::esgplotbands(sim_SVJD, main="SVJD Paths under the Physical Measure", 
                         xlab="Time",
                         ylab="Stock prices")
esgtoolkit::esgplotbands(sim_Heston, main="Heston Paths under the Physical Measure", 
                         xlab="Time",
                         ylab="Stock prices")                         

# Summary statistics
cat("Physical measure statistics (GBM):\n")
## Physical measure statistics (GBM):
terminal_prices_physical_GBM <- sim_GBM[nrow(sim_GBM), ]
cat("Mean terminal price:", mean(terminal_prices_physical_GBM), "\n")
## Mean terminal price: 150.3905
cat("Std terminal price:", sd(terminal_prices_physical_GBM), "\n")
## Std terminal price: 13.40266
cat("Expected under P:", S0 * exp(mu * maturity), "\n\n")
## Expected under P: 149.1825
cat("Physical measure statistics (SVJD):\n")
## Physical measure statistics (SVJD):
terminal_prices_physical_SVJD <- sim_SVJD[nrow(sim_SVJD), ]
cat("Mean terminal price:", mean(terminal_prices_physical_SVJD), "\n")
## Mean terminal price: 150.5358
cat("Std terminal price:", sd(terminal_prices_physical_SVJD), "\n")
## Std terminal price: 18.01574
cat("Expected under P:", S0 * exp(mu * maturity), "\n\n")
## Expected under P: 149.1825
cat("Physical measure statistics (Heston):\n")
## Physical measure statistics (Heston):
terminal_prices_physical_Heston <- sim_Heston[nrow(sim_Heston), ]
cat("Mean terminal price:", mean(terminal_prices_physical_Heston), "\n")
## Mean terminal price: 151.4444
cat("Std terminal price:", sd(terminal_prices_physical_Heston), "\n")
## Std terminal price: 16.65026
cat("Expected under P:", S0 * exp(mu * maturity), "\n\n")
## Expected under P: 149.1825
# =============================================================================
# 2 - COMPUTE DISCOUNTED PRICES (Transform to Martingale Domain)
# =============================================================================

discounted_prices_GBM <- esgtoolkit::esgdiscountfactor(r=r, X=sim_GBM)
discounted_prices_SVJD <- esgtoolkit::esgdiscountfactor(r=r, X=sim_SVJD)
discounted_prices_Heston <- esgtoolkit::esgdiscountfactor(r=r, X=sim_Heston)
martingale_diff_GBM <- discounted_prices_GBM - S0 
martingale_diff_SVJD <- discounted_prices_SVJD - S0 
martingale_diff_Heston <- discounted_prices_Heston - S0 

cat("Martingale differences dimensions (GBM):", dim(martingale_diff_GBM), "\n")
## Martingale differences dimensions (GBM): 1261 100
cat("Mean martingale diff (should be ≠ 0 under P):\n")
## Mean martingale diff (should be ≠ 0 under P):
print(t.test(rowMeans(martingale_diff_GBM)))
## 
##  One Sample t-test
## 
## data:  rowMeans(martingale_diff_GBM)
## t = 60.9, df = 1260, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  8.356654 8.912984
## sample estimates:
## mean of x 
##  8.634819
cat("\nMartingale differences dimensions (SVJD):", dim(martingale_diff_SVJD), "\n")
## 
## Martingale differences dimensions (SVJD): 1261 100
cat("Mean martingale diff (should be ≠ 0 under P):\n")
## Mean martingale diff (should be ≠ 0 under P):
print(t.test(rowMeans(martingale_diff_SVJD)))
## 
##  One Sample t-test
## 
## data:  rowMeans(martingale_diff_SVJD)
## t = 57.891, df = 1260, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  8.102235 8.670648
## sample estimates:
## mean of x 
##  8.386441
cat("\nMartingale differences dimensions (Heston):", dim(martingale_diff_Heston), "\n")
## 
## Martingale differences dimensions (Heston): 1261 100
cat("Mean martingale diff (should be ≠ 0 under P):\n")
## Mean martingale diff (should be ≠ 0 under P):
print(t.test(rowMeans(martingale_diff_Heston)))
## 
##  One Sample t-test
## 
## data:  rowMeans(martingale_diff_Heston)
## t = 58.357, df = 1260, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  8.384831 8.968209
## sample estimates:
## mean of x 
##   8.67652
# =============================================================================
# 3 - VISUALIZE RISK PREMIUM
# =============================================================================

par(mfrow=c(2,2))

matplot(discounted_prices_GBM, type='l', col=rgb(0,0,1,0.3),
        main="Discounted Stock Prices (Physical Measure - GBM)",
        ylab="exp(-rt) * S_t", xlab="Time step")
abline(h=S0, col='red', lwd=2, lty=2)

matplot(discounted_prices_SVJD, type='l', col=rgb(0,0,1,0.3),
        main="Discounted Stock Prices (Physical Measure - SVJD)",
        ylab="exp(-rt) * S_t", xlab="Time step")
abline(h=S0, col='red', lwd=2, lty=2)

matplot(discounted_prices_Heston, type='l', col=rgb(0,0,1,0.3),
        main="Discounted Stock Prices (Physical Measure - Heston)",
        ylab="exp(-rt) * S_t", xlab="Time step")
abline(h=S0, col='red', lwd=2, lty=2)

par(mfrow=c(1,1))

mean_disc_path_GBM <- rowMeans(discounted_prices_GBM)
times_plot <- as.numeric(time(discounted_prices_GBM)) 
plot(times_plot, mean_disc_path_GBM, type='l', lwd=2, col='blue',
     main="Risk Premium in Discounted Prices (GBM)",
     xlab="Time (years)", ylab="E[exp(-rt)*S_t]")
abline(h=S0, col='red', lwd=2, lty=2)
lines(times_plot, S0 * exp((mu-r)*times_plot), col='green', lwd=2, lty=3)
legend("topleft",
       legend=c("Empirical mean", "S0", "Theoretical (μ-r drift)"),
       col=c('blue','red','green'), lty=c(1,2,3), lwd=2)

mean_disc_path_SVJD <- rowMeans(discounted_prices_SVJD)
times_plot <- as.numeric(time(discounted_prices_SVJD)) 
plot(times_plot, mean_disc_path_SVJD, type='l', lwd=2, col='blue',
     main="Risk Premium in Discounted Prices (SVJD)",
     xlab="Time (years)", ylab="E[exp(-rt)*S_t]")
abline(h=S0, col='red', lwd=2, lty=2)
lines(times_plot, S0 * exp((mu-r)*times_plot), col='green', lwd=2, lty=3)
legend("topleft",
       legend=c("Empirical mean", "S0", "Theoretical (μ-r drift)"),
       col=c('blue','red','green'), lty=c(1,2,3), lwd=2)

mean_disc_path_Heston <- rowMeans(discounted_prices_Heston)
times_plot <- as.numeric(time(discounted_prices_Heston)) 
plot(times_plot, mean_disc_path_Heston, type='l', lwd=2, col='blue',
     main="Risk Premium in Discounted Prices (Heston)",
     xlab="Time (years)", ylab="E[exp(-rt)*S_t]")
abline(h=S0, col='red', lwd=2, lty=2)
lines(times_plot, S0 * exp((mu-r)*times_plot), col='green', lwd=2, lty=3)
legend("topleft",
       legend=c("Empirical mean", "S0", "Theoretical (μ-r drift)"),
       col=c('blue','red','green'), lty=c(1,2,3), lwd=2)     

# =============================================================================
# 4 - FIT ARIMA MODELS TO EXTRACT RISK PREMIUM
# =============================================================================

n_periods <- nrow(martingale_diff_GBM)
n_paths <- ncol(martingale_diff_GBM)

martingale_increments_GBM <- diff(martingale_diff_GBM)
martingale_increments_SVJD <- diff(martingale_diff_SVJD)
martingale_increments_Heston <- diff(martingale_diff_Heston)

# Initialize storage arrays
arima_residuals_GBM <- array(NA, dim = c(nrow(martingale_increments_GBM), n_paths))
centered_arima_residuals_GBM <- array(NA, dim = c(nrow(martingale_increments_GBM), n_paths))
means_arima_residuals_GBM <- rep(NA, n_paths)
arima_models_GBM <- list()

arima_residuals_SVJD <- array(NA, dim = c(nrow(martingale_increments_SVJD), n_paths))
centered_arima_residuals_SVJD <- array(NA, dim = c(nrow(martingale_increments_SVJD), n_paths))
means_arima_residuals_SVJD <- rep(NA, n_paths)
arima_models_SVJD <- list()

arima_residuals_Heston <- array(NA, dim = c(nrow(martingale_increments_Heston), n_paths))
centered_arima_residuals_Heston <- array(NA, dim = c(nrow(martingale_increments_Heston), n_paths))
means_arima_residuals_Heston <- rep(NA, n_paths)
arima_models_Heston <- list()

# Fit ARIMA models to GBM
cat("Fitting ARIMA models to", n_paths, "GBM paths...\n")
## Fitting ARIMA models to 100 GBM paths...
for (i in 1:n_paths) {
  y <- as.numeric(martingale_increments_GBM[, i])
  fit <- forecast::thetaf(y)
  arima_models_GBM[[i]] <- fit
  
  res <- as.numeric(residuals(fit))
  arima_residuals_GBM[, i] <- res
  
  centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
  means_arima_residuals_GBM[i] <- attr(centre_arima_residuals, "scaled:center")
  centered_arima_residuals_GBM[, i] <- centre_arima_residuals[,1]
}

# Fit ARIMA models to SVJD
cat("Fitting ARIMA models to", n_paths, "SVJD paths...\n")
## Fitting ARIMA models to 100 SVJD paths...
for (i in 1:n_paths) {
  y <- as.numeric(martingale_increments_SVJD[, i])
  fit <- forecast::thetaf(y)
  arima_models_SVJD[[i]] <- fit
  
  res <- as.numeric(residuals(fit))
  arima_residuals_SVJD[, i] <- res
  
  centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
  means_arima_residuals_SVJD[i] <- attr(centre_arima_residuals, "scaled:center")
  centered_arima_residuals_SVJD[, i] <- centre_arima_residuals[,1]
}

# Fit ARIMA models to Heston
cat("Fitting ARIMA models to", n_paths, "Heston paths...\n")
## Fitting ARIMA models to 100 Heston paths...
for (i in 1:n_paths) {
  y <- as.numeric(martingale_increments_Heston[, i])
  fit <- forecast::thetaf(y)
  arima_models_Heston[[i]] <- fit
  
  res <- as.numeric(residuals(fit))
  arima_residuals_Heston[, i] <- res
  
  centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
  means_arima_residuals_Heston[i] <- attr(centre_arima_residuals, "scaled:center")
  centered_arima_residuals_Heston[, i] <- centre_arima_residuals[,1]
}

# Box-Ljung tests
pvalues_GBM <- sapply(1:ncol(centered_arima_residuals_GBM), 
                  function(i) Box.test(centered_arima_residuals_GBM[,i])$p.value)
cat("\nBox-Ljung test p-values (GBM):\n")
## 
## Box-Ljung test p-values (GBM):
cat("Mean p-value:", mean(pvalues_GBM), "\n")
## Mean p-value: 0.5035183
cat("Proportion > 0.05:", mean(pvalues_GBM > 0.05), "\n")
## Proportion > 0.05: 0.96
pvalues_SVJD <- sapply(1:ncol(centered_arima_residuals_SVJD), 
                  function(i) Box.test(centered_arima_residuals_SVJD[,i])$p.value)
cat("\nBox-Ljung test p-values (SVJD):\n")
## 
## Box-Ljung test p-values (SVJD):
cat("Mean p-value:", mean(pvalues_SVJD), "\n")
## Mean p-value: 0.4609734
cat("Proportion > 0.05:", mean(pvalues_SVJD > 0.05), "\n")
## Proportion > 0.05: 0.84
pvalues_Heston <- sapply(1:ncol(centered_arima_residuals_Heston), 
                  function(i) Box.test(centered_arima_residuals_Heston[,i])$p.value)
cat("\nBox-Ljung test p-values (Heston):\n")
## 
## Box-Ljung test p-values (Heston):
cat("Mean p-value:", mean(pvalues_Heston), "\n")
## Mean p-value: 0.4114637
cat("Proportion > 0.05:", mean(pvalues_Heston > 0.05), "\n")
## Proportion > 0.05: 0.81
par(mfrow=c(1,3))
hist(pvalues_GBM, breaks=20, col='lightgreen',
     main="Box-Ljung P-values (GBM)",
     xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)

hist(pvalues_SVJD, breaks=20, col='lightblue',
     main="Box-Ljung P-values (SVJD)",
     xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)

hist(pvalues_Heston, breaks=20, col='lightcoral',
     main="Box-Ljung P-values (Heston)",
     xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)

par(mfrow=c(1,1))
# =============================================================================
# 5 - GENERATE RISK-NEUTRAL PATHS
# =============================================================================

bootstrap_independent <- function(x, R = 1000, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  
  n <- length(x)
  
  # Each column is an independent bootstrap sample
  replicate(R, sample(x, size = n, replace = TRUE))
}

cat("\n\nGenerating risk-neutral paths from ALL historical paths...\n")
## 
## 
## Generating risk-neutral paths from ALL historical paths...
n_sim_per_path <- 20  # Generate 10 paths per historical path
times <- seq(0, maturity, length.out = n_periods)
discount_factor <- exp(r * times)

# Storage for all risk-neutral paths
all_S_tilde_GBM <- list()
all_S_tilde_SVJD <- list()
all_S_tilde_Heston <- list()

# Generate GBM risk-neutral paths
for (i in 1:n_paths) {
  resampled_residuals <- bootstrap_independent(centered_arima_residuals_GBM[, i], 
                                               R = n_sim_per_path)
  fit <- arima_models_GBM[[i]]
  fitted_increments <- as.numeric(fitted(fit))
  
  discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
  discounted_path[1, ] <- S0
  
  increments <- matrix(scale(fitted_increments, center = TRUE, 
                             scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) + 
                resampled_residuals[1:(n_periods - 1), ]
  
  discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
  S_tilde_price <- discounted_path * discount_factor
  all_S_tilde_GBM[[i]] <- S_tilde_price
}

# Generate SVJD risk-neutral paths
for (i in 1:n_paths) {
  resampled_residuals <- bootstrap_independent(centered_arima_residuals_SVJD[, i], 
                                               R = n_sim_per_path)
  fit <- arima_models_SVJD[[i]]
  fitted_increments <- as.numeric(fitted(fit))
  
  discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
  discounted_path[1, ] <- S0
  
  increments <- matrix(scale(fitted_increments, center = TRUE, 
                             scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) + 
                resampled_residuals[1:(n_periods - 1), ]
  
  discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
  S_tilde_price <- discounted_path * discount_factor
  all_S_tilde_SVJD[[i]] <- S_tilde_price
}

# Generate Heston risk-neutral paths
for (i in 1:n_paths) {
  resampled_residuals <- bootstrap_independent(centered_arima_residuals_Heston[, i], 
                                               R = n_sim_per_path)
  fit <- arima_models_Heston[[i]]
  fitted_increments <- as.numeric(fitted(fit))
  
  discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
  discounted_path[1, ] <- S0
  
  increments <- matrix(scale(fitted_increments, center = TRUE, 
                             scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) + 
                resampled_residuals[1:(n_periods - 1), ]
  
  discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
  S_tilde_price <- discounted_path * discount_factor
  all_S_tilde_Heston[[i]] <- S_tilde_price
}

# Combine all paths
S_tilde_combined_GBM <- do.call(cbind, all_S_tilde_GBM)
cat("Total GBM risk-neutral paths generated:", ncol(S_tilde_combined_GBM), "\n")
## Total GBM risk-neutral paths generated: 2000
S_tilde_combined_SVJD <- do.call(cbind, all_S_tilde_SVJD)
cat("Total SVJD risk-neutral paths generated:", ncol(S_tilde_combined_SVJD), "\n")
## Total SVJD risk-neutral paths generated: 2000
S_tilde_combined_Heston <- do.call(cbind, all_S_tilde_Heston)
cat("Total Heston risk-neutral paths generated:", ncol(S_tilde_combined_Heston), "\n\n")
## Total Heston risk-neutral paths generated: 2000
# Convert to time series
S_tilde_ts_GBM <- ts(S_tilde_combined_GBM, start = start(sim_GBM), 
                     frequency = frequency(sim_GBM))
S_tilde_ts_SVJD <- ts(S_tilde_combined_SVJD, start = start(sim_SVJD), 
                      frequency = frequency(sim_SVJD))
S_tilde_ts_Heston <- ts(S_tilde_combined_Heston, start = start(sim_Heston), 
                        frequency = frequency(sim_Heston))

# Visualize risk-neutral paths
par(mfrow=c(1,3))
esgtoolkit::esgplotbands(S_tilde_ts_GBM, 
                        main="Risk-Neutral Paths - GBM")
esgtoolkit::esgplotbands(S_tilde_ts_SVJD, 
                        main="Risk-Neutral Paths - SVJD")
esgtoolkit::esgplotbands(S_tilde_ts_Heston, 
                        main="Risk-Neutral Paths - Heston")

# Sample plots
par(mfrow=c(1,3))
matplot(S_tilde_combined_GBM[, sample(ncol(S_tilde_combined_GBM), 200)], 
        type='l', col=rgb(0,0,1,0.1),
        main="Risk-Neutral Paths (GBM)",
        xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_GBM), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)

matplot(S_tilde_combined_SVJD[, sample(ncol(S_tilde_combined_SVJD), 200)], 
        type='l', col=rgb(0,0,1,0.1),
        main="Risk-Neutral Paths (SVJD)",
        xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_SVJD), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)

matplot(S_tilde_combined_Heston[, sample(ncol(S_tilde_combined_Heston), 200)], 
        type='l', col=rgb(0,0,1,0.1),
        main="Risk-Neutral Paths (Heston)",
        xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_Heston), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)

par(mfrow=c(1,1))
# =============================================================================
# 6 - VERIFY RISK-NEUTRAL PROPERTY
# =============================================================================

cat("\n=== RISK-NEUTRAL VERIFICATION ===\n\n")
## 
## === RISK-NEUTRAL VERIFICATION ===
terminal_prices_rn_GBM <- S_tilde_combined_GBM[n_periods, ]
terminal_prices_rn_SVJD <- S_tilde_combined_SVJD[n_periods, ]
terminal_prices_rn_Heston <- S_tilde_combined_Heston[n_periods, ]
capitalized_stock_price <- S0 * exp(r * maturity)

cat("GBM Risk-Neutral Verification:\n")
## GBM Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_GBM), "\n")
## Empirical mean: 128.3881
cat("Difference:", mean(terminal_prices_rn_GBM) - capitalized_stock_price, "\n")
## Difference: -0.01439695
print(t.test(terminal_prices_rn_GBM - capitalized_stock_price))
## 
##  One Sample t-test
## 
## data:  terminal_prices_rn_GBM - capitalized_stock_price
## t = -0.051539, df = 1999, p-value = 0.9589
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.5622233  0.5334294
## sample estimates:
##   mean of x 
## -0.01439695
cat("\nSVJD Risk-Neutral Verification:\n")
## 
## SVJD Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_SVJD), "\n")
## Empirical mean: 128.1574
cat("Difference:", mean(terminal_prices_rn_SVJD) - capitalized_stock_price, "\n")
## Difference: -0.2451161
print(t.test(terminal_prices_rn_SVJD - capitalized_stock_price))
## 
##  One Sample t-test
## 
## data:  terminal_prices_rn_SVJD - capitalized_stock_price
## t = -0.61866, df = 1999, p-value = 0.5362
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.0221290  0.5318968
## sample estimates:
##  mean of x 
## -0.2451161
cat("\nHeston Risk-Neutral Verification:\n")
## 
## Heston Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_Heston), "\n")
## Empirical mean: 128.5378
cat("Difference:", mean(terminal_prices_rn_Heston) - capitalized_stock_price, "\n")
## Difference: 0.1352577
print(t.test(terminal_prices_rn_Heston - capitalized_stock_price))
## 
##  One Sample t-test
## 
## data:  terminal_prices_rn_Heston - capitalized_stock_price
## t = 0.40195, df = 1999, p-value = 0.6878
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.5246852  0.7952007
## sample estimates:
## mean of x 
## 0.1352577
# Visualization comparison
par(mfrow=c(3, 2))
hist(terminal_prices_physical_GBM, breaks=30, col=rgb(1,0,0,0.5),
     main="Terminal Prices: Physical (GBM)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_GBM), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)

hist(terminal_prices_rn_GBM, breaks=30, col=rgb(0,0,1,0.5),
     main="Terminal Prices: Risk-Neutral (GBM)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_GBM), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)

hist(terminal_prices_physical_SVJD, breaks=30, col=rgb(1,0,0,0.5),
     main="Terminal Prices: Physical (SVJD)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_SVJD), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)

hist(terminal_prices_rn_SVJD, breaks=30, col=rgb(0,0,1,0.5),
     main="Terminal Prices: Risk-Neutral (SVJD)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_SVJD), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)

hist(terminal_prices_physical_Heston, breaks=30, col=rgb(1,0,0,0.5),
     main="Terminal Prices: Physical (Heston)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_Heston), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)

hist(terminal_prices_rn_Heston, breaks=30, col=rgb(0,0,1,0.5),
     main="Terminal Prices: Risk-Neutral (Heston)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_Heston), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)

par(mfrow=c(1,1))
# =============================================================================
# 7 - OPTION PRICING
# =============================================================================

cat("\n=== OPTION PRICING ===\n\n")
## 
## === OPTION PRICING ===
bs_price <- function(S, K, r, sigma, T, q = 0) {

  d1 <- (log(S / K) + (r - q + 0.5 * sigma^2) * T) / (sigma * sqrt(T))
  d2 <- d1 - sigma * sqrt(T)

  call <- S * exp(-q * T) * pnorm(d1) - K * exp(-r * T) * pnorm(d2)
  put  <- K * exp(-r * T) * pnorm(-d2) - S * exp(-q * T) * pnorm(-d1)

  list(call = call, put = put)
}


strikes <- seq(80, 160, by=10)
d_f <- exp(-r * maturity)

# Function to price options
price_options <- function(terminal_prices, strikes, discount_factor) {
  n_strikes <- length(strikes)
  call_prices <- numeric(n_strikes)
  bs_call_prices <- numeric(n_strikes)
  put_prices <- numeric(n_strikes)
  bs_put_prices <- numeric(n_strikes)
  call_se <- numeric(n_strikes)
  put_se <- numeric(n_strikes)
  
  for (k in 1:n_strikes) {
    K <- strikes[k]
    call_payoffs <- pmax(terminal_prices - K, 0)
    call_prices[k] <- mean(call_payoffs) * discount_factor
    call_se[k] <- sd(call_payoffs) / sqrt(length(call_payoffs)) * discount_factor
    
    put_payoffs <- pmax(K - terminal_prices, 0)
    put_prices[k] <- mean(put_payoffs) * discount_factor
    put_se[k] <- sd(put_payoffs) / sqrt(length(put_payoffs)) * discount_factor
    
    bs_prices <- bs_price(S0, K, r, sigma=sigma, T=5, q = 0)
    bs_call_prices[k] <- bs_prices$call
    bs_put_prices[k] <- bs_prices$put
  }
  
  list(call_prices = call_prices, put_prices = put_prices,
       bs_call_prices = bs_call_prices, bs_put_prices = bs_put_prices,
       call_se = call_se, put_se = put_se)
}

# Price options for all models
options_GBM <- price_options(terminal_prices_rn_GBM, strikes, d_f)
options_SVJD <- price_options(terminal_prices_rn_SVJD, strikes, d_f)
options_Heston <- price_options(terminal_prices_rn_Heston, strikes, d_f)
kableExtra::kable(as.data.frame(options_GBM))
call_prices put_prices bs_call_prices bs_put_prices call_se put_se
37.6847250 0.0000000 37.6959374 0.0000001 0.2175495 0.0000000
29.8980644 0.0013472 29.9079897 0.0000602 0.2174527 0.0013472
22.1452991 0.0365898 22.1260239 0.0061022 0.2154791 0.0092665
14.6435444 0.3228429 14.4726589 0.1407450 0.2036358 0.0336206
8.0122810 1.4795873 7.6644048 1.1204988 0.1710185 0.0788649
3.2699519 4.5252661 3.0014135 4.2455153 0.1165797 0.1375942
0.9066440 9.9499659 0.8280730 9.8601826 0.0606846 0.1860634
0.1685717 16.9999015 0.1608817 16.9809992 0.0227524 0.2096258
0.0121927 24.6315303 0.0226357 24.6307609 0.0053122 0.2167928
kableExtra::kable(as.data.frame(options_Heston))
call_prices put_prices bs_call_prices bs_put_prices call_se put_se
37.8028720 0.0015958 37.6959374 0.0000001 0.2619526 0.0015958
30.0399788 0.0267105 29.9079897 0.0000602 0.2603824 0.0089609
22.3576823 0.1324218 22.1260239 0.0061022 0.2552281 0.0240528
15.0348491 0.5975964 14.4726589 0.1407450 0.2383996 0.0534645
8.6676523 2.0184075 7.6644048 1.1204988 0.2021327 0.1015914
4.0766962 5.2154591 3.0014135 4.2455153 0.1478793 0.1598164
1.5391986 10.4659694 0.8280730 9.8601826 0.0911697 0.2103637
0.4459361 17.1607147 0.1608817 16.9809992 0.0486331 0.2421993
0.1148506 24.6176371 0.0226357 24.6307609 0.0241640 0.2554788
kableExtra::kable(as.data.frame(options_SVJD))
call_prices put_prices bs_call_prices bs_put_prices call_se put_se
37.5880705 0.0830298 37.6959374 0.0000001 0.3025649 0.0233010
29.9091959 0.1921630 29.9079897 0.0000602 0.2964688 0.0395838
22.3743079 0.4452828 22.1260239 0.0061022 0.2853585 0.0617545
15.2673169 1.1262996 14.4726589 0.1407450 0.2627206 0.0947876
9.1376458 2.7846363 7.6644048 1.1204988 0.2229804 0.1415378
4.5976251 6.0326235 3.0014135 4.2455153 0.1689933 0.1972375
1.9332399 11.1562461 0.8280730 9.8601826 0.1140386 0.2462269
0.7097549 17.7207690 0.1608817 16.9809992 0.0705917 0.2786469
0.2313729 25.0303948 0.0226357 24.6307609 0.0434949 0.2958461

ETS

# =============================================================================
# 4 - FIT ARIMA MODELS TO EXTRACT RISK PREMIUM
# =============================================================================

n_periods <- nrow(martingale_diff_GBM)
n_paths <- ncol(martingale_diff_GBM)

martingale_increments_GBM <- diff(martingale_diff_GBM)
martingale_increments_SVJD <- diff(martingale_diff_SVJD)
martingale_increments_Heston <- diff(martingale_diff_Heston)

# Initialize storage arrays
arima_residuals_GBM <- array(NA, dim = c(nrow(martingale_increments_GBM), n_paths))
centered_arima_residuals_GBM <- array(NA, dim = c(nrow(martingale_increments_GBM), n_paths))
means_arima_residuals_GBM <- rep(NA, n_paths)
arima_models_GBM <- list()

arima_residuals_SVJD <- array(NA, dim = c(nrow(martingale_increments_SVJD), n_paths))
centered_arima_residuals_SVJD <- array(NA, dim = c(nrow(martingale_increments_SVJD), n_paths))
means_arima_residuals_SVJD <- rep(NA, n_paths)
arima_models_SVJD <- list()

arima_residuals_Heston <- array(NA, dim = c(nrow(martingale_increments_Heston), n_paths))
centered_arima_residuals_Heston <- array(NA, dim = c(nrow(martingale_increments_Heston), n_paths))
means_arima_residuals_Heston <- rep(NA, n_paths)
arima_models_Heston <- list()

# Fit ARIMA models to GBM
cat("Fitting ARIMA models to", n_paths, "GBM paths...\n")
## Fitting ARIMA models to 100 GBM paths...
for (i in 1:n_paths) {
  y <- as.numeric(martingale_increments_GBM[, i])
  fit <- forecast::ets(y)
  arima_models_GBM[[i]] <- fit
  
  res <- as.numeric(residuals(fit))
  arima_residuals_GBM[, i] <- res
  
  centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
  means_arima_residuals_GBM[i] <- attr(centre_arima_residuals, "scaled:center")
  centered_arima_residuals_GBM[, i] <- centre_arima_residuals[,1]
}

# Fit ARIMA models to SVJD
cat("Fitting ARIMA models to", n_paths, "SVJD paths...\n")
## Fitting ARIMA models to 100 SVJD paths...
for (i in 1:n_paths) {
  y <- as.numeric(martingale_increments_SVJD[, i])
  fit <- forecast::ets(y)
  arima_models_SVJD[[i]] <- fit
  
  res <- as.numeric(residuals(fit))
  arima_residuals_SVJD[, i] <- res
  
  centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
  means_arima_residuals_SVJD[i] <- attr(centre_arima_residuals, "scaled:center")
  centered_arima_residuals_SVJD[, i] <- centre_arima_residuals[,1]
}

# Fit ARIMA models to Heston
cat("Fitting ARIMA models to", n_paths, "Heston paths...\n")
## Fitting ARIMA models to 100 Heston paths...
for (i in 1:n_paths) {
  y <- as.numeric(martingale_increments_Heston[, i])
  fit <- forecast::thetaf(y)
  arima_models_Heston[[i]] <- fit
  
  res <- as.numeric(residuals(fit))
  arima_residuals_Heston[, i] <- res
  
  centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
  means_arima_residuals_Heston[i] <- attr(centre_arima_residuals, "scaled:center")
  centered_arima_residuals_Heston[, i] <- centre_arima_residuals[,1]
}

# Box-Ljung tests
pvalues_GBM <- sapply(1:ncol(centered_arima_residuals_GBM), 
                  function(i) Box.test(centered_arima_residuals_GBM[,i])$p.value)
cat("\nBox-Ljung test p-values (GBM):\n")
## 
## Box-Ljung test p-values (GBM):
cat("Mean p-value:", mean(pvalues_GBM), "\n")
## Mean p-value: 0.5048272
cat("Proportion > 0.05:", mean(pvalues_GBM > 0.05), "\n")
## Proportion > 0.05: 0.96
pvalues_SVJD <- sapply(1:ncol(centered_arima_residuals_SVJD), 
                  function(i) Box.test(centered_arima_residuals_SVJD[,i])$p.value)
cat("\nBox-Ljung test p-values (SVJD):\n")
## 
## Box-Ljung test p-values (SVJD):
cat("Mean p-value:", mean(pvalues_SVJD), "\n")
## Mean p-value: 0.4611386
cat("Proportion > 0.05:", mean(pvalues_SVJD > 0.05), "\n")
## Proportion > 0.05: 0.84
pvalues_Heston <- sapply(1:ncol(centered_arima_residuals_Heston), 
                  function(i) Box.test(centered_arima_residuals_Heston[,i])$p.value)
cat("\nBox-Ljung test p-values (Heston):\n")
## 
## Box-Ljung test p-values (Heston):
cat("Mean p-value:", mean(pvalues_Heston), "\n")
## Mean p-value: 0.4114637
cat("Proportion > 0.05:", mean(pvalues_Heston > 0.05), "\n")
## Proportion > 0.05: 0.81
par(mfrow=c(1,3))
hist(pvalues_GBM, breaks=20, col='lightgreen',
     main="Box-Ljung P-values (GBM)",
     xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)

hist(pvalues_SVJD, breaks=20, col='lightblue',
     main="Box-Ljung P-values (SVJD)",
     xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)

hist(pvalues_Heston, breaks=20, col='lightcoral',
     main="Box-Ljung P-values (Heston)",
     xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)

par(mfrow=c(1,1))
# =============================================================================
# 5 - GENERATE RISK-NEUTRAL PATHS
# =============================================================================

bootstrap_independent <- function(x, R = 1000, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  
  n <- length(x)
  
  # Each column is an independent bootstrap sample
  replicate(R, sample(x, size = n, replace = TRUE))
}

cat("\n\nGenerating risk-neutral paths from ALL historical paths...\n")
## 
## 
## Generating risk-neutral paths from ALL historical paths...
n_sim_per_path <- 20  # Generate 10 paths per historical path
times <- seq(0, maturity, length.out = n_periods)
discount_factor <- exp(r * times)

# Storage for all risk-neutral paths
all_S_tilde_GBM <- list()
all_S_tilde_SVJD <- list()
all_S_tilde_Heston <- list()

# Generate GBM risk-neutral paths
for (i in 1:n_paths) {
  resampled_residuals <- bootstrap_independent(centered_arima_residuals_GBM[, i], 
                                               R = n_sim_per_path)
  fit <- arima_models_GBM[[i]]
  fitted_increments <- as.numeric(fitted(fit))
  
  discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
  discounted_path[1, ] <- S0
  
  increments <- matrix(scale(fitted_increments, center = TRUE, 
                             scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) + 
                resampled_residuals[1:(n_periods - 1), ]
  
  discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
  S_tilde_price <- discounted_path * discount_factor
  all_S_tilde_GBM[[i]] <- S_tilde_price
}

# Generate SVJD risk-neutral paths
for (i in 1:n_paths) {
  resampled_residuals <- bootstrap_independent(centered_arima_residuals_SVJD[, i], 
                                               R = n_sim_per_path)
  fit <- arima_models_SVJD[[i]]
  fitted_increments <- as.numeric(fitted(fit))
  
  discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
  discounted_path[1, ] <- S0
  
  increments <- matrix(scale(fitted_increments, center = TRUE, 
                             scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) + 
                resampled_residuals[1:(n_periods - 1), ]
  
  discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
  S_tilde_price <- discounted_path * discount_factor
  all_S_tilde_SVJD[[i]] <- S_tilde_price
}

# Generate Heston risk-neutral paths
for (i in 1:n_paths) {
  resampled_residuals <- bootstrap_independent(centered_arima_residuals_Heston[, i], 
                                               R = n_sim_per_path)
  fit <- arima_models_Heston[[i]]
  fitted_increments <- as.numeric(fitted(fit))
  
  discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
  discounted_path[1, ] <- S0
  
  increments <- matrix(scale(fitted_increments, center = TRUE, 
                             scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) + 
                resampled_residuals[1:(n_periods - 1), ]
  
  discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
  S_tilde_price <- discounted_path * discount_factor
  all_S_tilde_Heston[[i]] <- S_tilde_price
}

# Combine all paths
S_tilde_combined_GBM <- do.call(cbind, all_S_tilde_GBM)
cat("Total GBM risk-neutral paths generated:", ncol(S_tilde_combined_GBM), "\n")
## Total GBM risk-neutral paths generated: 2000
S_tilde_combined_SVJD <- do.call(cbind, all_S_tilde_SVJD)
cat("Total SVJD risk-neutral paths generated:", ncol(S_tilde_combined_SVJD), "\n")
## Total SVJD risk-neutral paths generated: 2000
S_tilde_combined_Heston <- do.call(cbind, all_S_tilde_Heston)
cat("Total Heston risk-neutral paths generated:", ncol(S_tilde_combined_Heston), "\n\n")
## Total Heston risk-neutral paths generated: 2000
# Convert to time series
S_tilde_ts_GBM <- ts(S_tilde_combined_GBM, start = start(sim_GBM), 
                     frequency = frequency(sim_GBM))
S_tilde_ts_SVJD <- ts(S_tilde_combined_SVJD, start = start(sim_SVJD), 
                      frequency = frequency(sim_SVJD))
S_tilde_ts_Heston <- ts(S_tilde_combined_Heston, start = start(sim_Heston), 
                        frequency = frequency(sim_Heston))

# Visualize risk-neutral paths
par(mfrow=c(1,3))
esgtoolkit::esgplotbands(S_tilde_ts_GBM, 
                        main="Risk-Neutral Paths - GBM")
esgtoolkit::esgplotbands(S_tilde_ts_SVJD, 
                        main="Risk-Neutral Paths - SVJD")
esgtoolkit::esgplotbands(S_tilde_ts_Heston, 
                        main="Risk-Neutral Paths - Heston")

# Sample plots
par(mfrow=c(1,3))
matplot(S_tilde_combined_GBM[, sample(ncol(S_tilde_combined_GBM), 200)], 
        type='l', col=rgb(0,0,1,0.1),
        main="Risk-Neutral Paths (GBM)",
        xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_GBM), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)

matplot(S_tilde_combined_SVJD[, sample(ncol(S_tilde_combined_SVJD), 200)], 
        type='l', col=rgb(0,0,1,0.1),
        main="Risk-Neutral Paths (SVJD)",
        xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_SVJD), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)

matplot(S_tilde_combined_Heston[, sample(ncol(S_tilde_combined_Heston), 200)], 
        type='l', col=rgb(0,0,1,0.1),
        main="Risk-Neutral Paths (Heston)",
        xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_Heston), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)

par(mfrow=c(1,1))
# =============================================================================
# 6 - VERIFY RISK-NEUTRAL PROPERTY
# =============================================================================

cat("\n=== RISK-NEUTRAL VERIFICATION ===\n\n")
## 
## === RISK-NEUTRAL VERIFICATION ===
terminal_prices_rn_GBM <- S_tilde_combined_GBM[n_periods, ]
terminal_prices_rn_SVJD <- S_tilde_combined_SVJD[n_periods, ]
terminal_prices_rn_Heston <- S_tilde_combined_Heston[n_periods, ]
capitalized_stock_price <- S0 * exp(r * maturity)

cat("GBM Risk-Neutral Verification:\n")
## GBM Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_GBM), "\n")
## Empirical mean: 127.924
cat("Difference:", mean(terminal_prices_rn_GBM) - capitalized_stock_price, "\n")
## Difference: -0.4785094
print(t.test(terminal_prices_rn_GBM - capitalized_stock_price))
## 
##  One Sample t-test
## 
## data:  terminal_prices_rn_GBM - capitalized_stock_price
## t = -1.7301, df = 1999, p-value = 0.08377
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.0209334  0.0639145
## sample estimates:
##  mean of x 
## -0.4785094
cat("\nSVJD Risk-Neutral Verification:\n")
## 
## SVJD Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_SVJD), "\n")
## Empirical mean: 127.8367
cat("Difference:", mean(terminal_prices_rn_SVJD) - capitalized_stock_price, "\n")
## Difference: -0.5658013
print(t.test(terminal_prices_rn_SVJD - capitalized_stock_price))
## 
##  One Sample t-test
## 
## data:  terminal_prices_rn_SVJD - capitalized_stock_price
## t = -1.3747, df = 1999, p-value = 0.1694
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.3729490  0.2413465
## sample estimates:
##  mean of x 
## -0.5658013
cat("\nHeston Risk-Neutral Verification:\n")
## 
## Heston Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_Heston), "\n")
## Empirical mean: 128.3782
cat("Difference:", mean(terminal_prices_rn_Heston) - capitalized_stock_price, "\n")
## Difference: -0.02433147
print(t.test(terminal_prices_rn_Heston - capitalized_stock_price))
## 
##  One Sample t-test
## 
## data:  terminal_prices_rn_Heston - capitalized_stock_price
## t = -0.072939, df = 1999, p-value = 0.9419
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.6785473  0.6298844
## sample estimates:
##   mean of x 
## -0.02433147
# Visualization comparison
par(mfrow=c(3, 2))
hist(terminal_prices_physical_GBM, breaks=30, col=rgb(1,0,0,0.5),
     main="Terminal Prices: Physical (GBM)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_GBM), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)

hist(terminal_prices_rn_GBM, breaks=30, col=rgb(0,0,1,0.5),
     main="Terminal Prices: Risk-Neutral (GBM)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_GBM), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)

hist(terminal_prices_physical_SVJD, breaks=30, col=rgb(1,0,0,0.5),
     main="Terminal Prices: Physical (SVJD)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_SVJD), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)

hist(terminal_prices_rn_SVJD, breaks=30, col=rgb(0,0,1,0.5),
     main="Terminal Prices: Risk-Neutral (SVJD)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_SVJD), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)

hist(terminal_prices_physical_Heston, breaks=30, col=rgb(1,0,0,0.5),
     main="Terminal Prices: Physical (Heston)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_Heston), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)

hist(terminal_prices_rn_Heston, breaks=30, col=rgb(0,0,1,0.5),
     main="Terminal Prices: Risk-Neutral (Heston)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_Heston), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)

par(mfrow=c(1,1))
# =============================================================================
# 7 - OPTION PRICING
# =============================================================================

cat("\n=== OPTION PRICING ===\n\n")
## 
## === OPTION PRICING ===
bs_price <- function(S, K, r, sigma, T, q = 0) {

  d1 <- (log(S / K) + (r - q + 0.5 * sigma^2) * T) / (sigma * sqrt(T))
  d2 <- d1 - sigma * sqrt(T)

  call <- S * exp(-q * T) * pnorm(d1) - K * exp(-r * T) * pnorm(d2)
  put  <- K * exp(-r * T) * pnorm(-d2) - S * exp(-q * T) * pnorm(-d1)

  list(call = call, put = put)
}


strikes <- seq(80, 160, by=10)
d_f <- exp(-r * maturity)

# Function to price options
price_options <- function(terminal_prices, strikes, discount_factor) {
  n_strikes <- length(strikes)
  call_prices <- numeric(n_strikes)
  bs_call_prices <- numeric(n_strikes)
  put_prices <- numeric(n_strikes)
  bs_put_prices <- numeric(n_strikes)
  call_se <- numeric(n_strikes)
  put_se <- numeric(n_strikes)
  
  for (k in 1:n_strikes) {
    K <- strikes[k]
    call_payoffs <- pmax(terminal_prices - K, 0)
    call_prices[k] <- mean(call_payoffs) * discount_factor
    call_se[k] <- sd(call_payoffs) / sqrt(length(call_payoffs)) * discount_factor
    
    put_payoffs <- pmax(K - terminal_prices, 0)
    put_prices[k] <- mean(put_payoffs) * discount_factor
    put_se[k] <- sd(put_payoffs) / sqrt(length(put_payoffs)) * discount_factor
    
    bs_prices <- bs_price(S0, K, r, sigma=sigma, T=5, q = 0)
    bs_call_prices[k] <- bs_prices$call
    bs_put_prices[k] <- bs_prices$put
  }
  
  list(call_prices = call_prices, put_prices = put_prices,
       bs_call_prices = bs_call_prices, bs_put_prices = bs_put_prices,
       call_se = call_se, put_se = put_se)
}

# Price options for all models
options_GBM <- price_options(terminal_prices_rn_GBM, strikes, d_f)
options_SVJD <- price_options(terminal_prices_rn_SVJD, strikes, d_f)
options_Heston <- price_options(terminal_prices_rn_Heston, strikes, d_f)
kableExtra::kable(as.data.frame(options_GBM))
call_prices put_prices bs_call_prices bs_put_prices call_se put_se
37.3233395 0.0000657 37.6959374 0.0000001 0.2153985 0.0000657
29.5427541 0.0074881 29.9079897 0.0000602 0.2148244 0.0053024
21.7860973 0.0388391 22.1260239 0.0061022 0.2130421 0.0128624
14.2477281 0.2884778 14.4726589 0.1407450 0.2029315 0.0332499
7.6572043 1.4859618 7.6644048 1.1204988 0.1702183 0.0777218
3.0857449 4.7025103 3.0014135 4.2455153 0.1154527 0.1362044
0.8569653 10.2617385 0.8280730 9.8601826 0.0608434 0.1841160
0.1686439 17.3614249 0.1608817 16.9809992 0.0251250 0.2069742
0.0200318 25.0008207 0.0226357 24.6307609 0.0091880 0.2140408
kableExtra::kable(as.data.frame(options_Heston))
call_prices put_prices bs_call_prices bs_put_prices call_se put_se
37.6822414 0.0052534 37.6959374 0.0000001 0.2593967 0.0032356
29.9177530 0.0287729 29.9079897 0.0000602 0.2578997 0.0110308
22.2663566 0.1653843 22.1260239 0.0061022 0.2510846 0.0277011
14.9579035 0.6449390 14.4726589 0.1407450 0.2334187 0.0579578
8.5691693 2.0442127 7.6644048 1.1204988 0.1971864 0.1052936
3.9707040 5.2337551 3.0014135 4.2455153 0.1421274 0.1627970
1.4344219 10.4854809 0.8280730 9.8601826 0.0850850 0.2126207
0.3989463 17.2380131 0.1608817 16.9809992 0.0416221 0.2426567
0.0714887 24.6985633 0.0226357 24.6307609 0.0170685 0.2558071
kableExtra::kable(as.data.frame(options_SVJD))
call_prices put_prices bs_call_prices bs_put_prices call_se put_se
37.3936091 0.1383182 37.6959374 0.0000001 0.3098695 0.0393111
29.7336475 0.2663644 29.9079897 0.0000602 0.3029494 0.0551093
22.2189875 0.5397123 22.1260239 0.0061022 0.2912961 0.0767336
15.1512978 1.2600304 14.4726589 0.1407450 0.2679371 0.1088503
9.0474757 2.9442161 7.6644048 1.1204988 0.2286590 0.1542833
4.5472247 6.2319730 3.0014135 4.2455153 0.1762522 0.2081396
1.9296726 11.4024287 0.8280730 9.8601826 0.1245120 0.2553862
0.7570340 18.0177979 0.1608817 16.9809992 0.0854171 0.2860005
0.3383181 25.3870899 0.0226357 24.6307609 0.0595079 0.3010065

ARIMA

# =============================================================================
# 4 - FIT ARIMA MODELS TO EXTRACT RISK PREMIUM
# =============================================================================

n_periods <- nrow(martingale_diff_GBM)
n_paths <- ncol(martingale_diff_GBM)

martingale_increments_GBM <- diff(martingale_diff_GBM)
martingale_increments_SVJD <- diff(martingale_diff_SVJD)
martingale_increments_Heston <- diff(martingale_diff_Heston)

# Initialize storage arrays
arima_residuals_GBM <- array(NA, dim = c(nrow(martingale_increments_GBM), n_paths))
centered_arima_residuals_GBM <- array(NA, dim = c(nrow(martingale_increments_GBM), n_paths))
means_arima_residuals_GBM <- rep(NA, n_paths)
arima_models_GBM <- list()

arima_residuals_SVJD <- array(NA, dim = c(nrow(martingale_increments_SVJD), n_paths))
centered_arima_residuals_SVJD <- array(NA, dim = c(nrow(martingale_increments_SVJD), n_paths))
means_arima_residuals_SVJD <- rep(NA, n_paths)
arima_models_SVJD <- list()

arima_residuals_Heston <- array(NA, dim = c(nrow(martingale_increments_Heston), n_paths))
centered_arima_residuals_Heston <- array(NA, dim = c(nrow(martingale_increments_Heston), n_paths))
means_arima_residuals_Heston <- rep(NA, n_paths)
arima_models_Heston <- list()

# Fit ARIMA models to GBM
cat("Fitting ARIMA models to", n_paths, "GBM paths...\n")
## Fitting ARIMA models to 100 GBM paths...
for (i in 1:n_paths) {
  y <- as.numeric(martingale_increments_GBM[, i])
  fit <- forecast::auto.arima(y)
  arima_models_GBM[[i]] <- fit
  
  res <- as.numeric(residuals(fit))
  arima_residuals_GBM[, i] <- res
  
  centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
  means_arima_residuals_GBM[i] <- attr(centre_arima_residuals, "scaled:center")
  centered_arima_residuals_GBM[, i] <- centre_arima_residuals[,1]
}

# Fit ARIMA models to SVJD
cat("Fitting ARIMA models to", n_paths, "SVJD paths...\n")
## Fitting ARIMA models to 100 SVJD paths...
for (i in 1:n_paths) {
  y <- as.numeric(martingale_increments_SVJD[, i])
  fit <- forecast::auto.arima(y)
  arima_models_SVJD[[i]] <- fit
  
  res <- as.numeric(residuals(fit))
  arima_residuals_SVJD[, i] <- res
  
  centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
  means_arima_residuals_SVJD[i] <- attr(centre_arima_residuals, "scaled:center")
  centered_arima_residuals_SVJD[, i] <- centre_arima_residuals[,1]
}

# Fit ARIMA models to Heston
cat("Fitting ARIMA models to", n_paths, "Heston paths...\n")
## Fitting ARIMA models to 100 Heston paths...
for (i in 1:n_paths) {
  y <- as.numeric(martingale_increments_Heston[, i])
  fit <- forecast::thetaf(y)
  arima_models_Heston[[i]] <- fit
  
  res <- as.numeric(residuals(fit))
  arima_residuals_Heston[, i] <- res
  
  centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
  means_arima_residuals_Heston[i] <- attr(centre_arima_residuals, "scaled:center")
  centered_arima_residuals_Heston[, i] <- centre_arima_residuals[,1]
}

# Box-Ljung tests
pvalues_GBM <- sapply(1:ncol(centered_arima_residuals_GBM), 
                  function(i) Box.test(centered_arima_residuals_GBM[,i])$p.value)
cat("\nBox-Ljung test p-values (GBM):\n")
## 
## Box-Ljung test p-values (GBM):
cat("Mean p-value:", mean(pvalues_GBM), "\n")
## Mean p-value: 0.6636805
cat("Proportion > 0.05:", mean(pvalues_GBM > 0.05), "\n")
## Proportion > 0.05: 0.95
pvalues_SVJD <- sapply(1:ncol(centered_arima_residuals_SVJD), 
                  function(i) Box.test(centered_arima_residuals_SVJD[,i])$p.value)
cat("\nBox-Ljung test p-values (SVJD):\n")
## 
## Box-Ljung test p-values (SVJD):
cat("Mean p-value:", mean(pvalues_SVJD), "\n")
## Mean p-value: 0.6953426
cat("Proportion > 0.05:", mean(pvalues_SVJD > 0.05), "\n")
## Proportion > 0.05: 0.98
pvalues_Heston <- sapply(1:ncol(centered_arima_residuals_Heston), 
                  function(i) Box.test(centered_arima_residuals_Heston[,i])$p.value)
cat("\nBox-Ljung test p-values (Heston):\n")
## 
## Box-Ljung test p-values (Heston):
cat("Mean p-value:", mean(pvalues_Heston), "\n")
## Mean p-value: 0.4114637
cat("Proportion > 0.05:", mean(pvalues_Heston > 0.05), "\n")
## Proportion > 0.05: 0.81
par(mfrow=c(1,3))
hist(pvalues_GBM, breaks=20, col='lightgreen',
     main="Box-Ljung P-values (GBM)",
     xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)

hist(pvalues_SVJD, breaks=20, col='lightblue',
     main="Box-Ljung P-values (SVJD)",
     xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)

hist(pvalues_Heston, breaks=20, col='lightcoral',
     main="Box-Ljung P-values (Heston)",
     xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)

par(mfrow=c(1,1))
# =============================================================================
# 5 - GENERATE RISK-NEUTRAL PATHS
# =============================================================================

bootstrap_independent <- function(x, R = 1000, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)
  
  n <- length(x)
  
  # Each column is an independent bootstrap sample
  replicate(R, sample(x, size = n, replace = TRUE))
}

cat("\n\nGenerating risk-neutral paths from ALL historical paths...\n")
## 
## 
## Generating risk-neutral paths from ALL historical paths...
n_sim_per_path <- 20  # Generate 10 paths per historical path
times <- seq(0, maturity, length.out = n_periods)
discount_factor <- exp(r * times)

# Storage for all risk-neutral paths
all_S_tilde_GBM <- list()
all_S_tilde_SVJD <- list()
all_S_tilde_Heston <- list()

# Generate GBM risk-neutral paths
for (i in 1:n_paths) {
  resampled_residuals <- bootstrap_independent(centered_arima_residuals_GBM[, i], 
                                               R = n_sim_per_path)
  fit <- arima_models_GBM[[i]]
  fitted_increments <- as.numeric(fitted(fit))
  
  discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
  discounted_path[1, ] <- S0
  
  increments <- matrix(scale(fitted_increments, center = TRUE, 
                             scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) + 
                resampled_residuals[1:(n_periods - 1), ]
  
  discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
  S_tilde_price <- discounted_path * discount_factor
  all_S_tilde_GBM[[i]] <- S_tilde_price
}

# Generate SVJD risk-neutral paths
for (i in 1:n_paths) {
  resampled_residuals <- bootstrap_independent(centered_arima_residuals_SVJD[, i], 
                                               R = n_sim_per_path)
  fit <- arima_models_SVJD[[i]]
  fitted_increments <- as.numeric(fitted(fit))
  
  discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
  discounted_path[1, ] <- S0
  
  increments <- matrix(scale(fitted_increments, center = TRUE, 
                             scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) + 
                resampled_residuals[1:(n_periods - 1), ]
  
  discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
  S_tilde_price <- discounted_path * discount_factor
  all_S_tilde_SVJD[[i]] <- S_tilde_price
}

# Generate Heston risk-neutral paths
for (i in 1:n_paths) {
  resampled_residuals <- bootstrap_independent(centered_arima_residuals_Heston[, i], 
                                               R = n_sim_per_path)
  fit <- arima_models_Heston[[i]]
  fitted_increments <- as.numeric(fitted(fit))
  
  discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
  discounted_path[1, ] <- S0
  
  increments <- matrix(scale(fitted_increments, center = TRUE, 
                             scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) + 
                resampled_residuals[1:(n_periods - 1), ]
  
  discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
  S_tilde_price <- discounted_path * discount_factor
  all_S_tilde_Heston[[i]] <- S_tilde_price
}

# Combine all paths
S_tilde_combined_GBM <- do.call(cbind, all_S_tilde_GBM)
cat("Total GBM risk-neutral paths generated:", ncol(S_tilde_combined_GBM), "\n")
## Total GBM risk-neutral paths generated: 2000
S_tilde_combined_SVJD <- do.call(cbind, all_S_tilde_SVJD)
cat("Total SVJD risk-neutral paths generated:", ncol(S_tilde_combined_SVJD), "\n")
## Total SVJD risk-neutral paths generated: 2000
S_tilde_combined_Heston <- do.call(cbind, all_S_tilde_Heston)
cat("Total Heston risk-neutral paths generated:", ncol(S_tilde_combined_Heston), "\n\n")
## Total Heston risk-neutral paths generated: 2000
# Convert to time series
S_tilde_ts_GBM <- ts(S_tilde_combined_GBM, start = start(sim_GBM), 
                     frequency = frequency(sim_GBM))
S_tilde_ts_SVJD <- ts(S_tilde_combined_SVJD, start = start(sim_SVJD), 
                      frequency = frequency(sim_SVJD))
S_tilde_ts_Heston <- ts(S_tilde_combined_Heston, start = start(sim_Heston), 
                        frequency = frequency(sim_Heston))

# Visualize risk-neutral paths
par(mfrow=c(1,3))
esgtoolkit::esgplotbands(S_tilde_ts_GBM, 
                        main="Risk-Neutral Paths - GBM")
esgtoolkit::esgplotbands(S_tilde_ts_SVJD, 
                        main="Risk-Neutral Paths - SVJD")
esgtoolkit::esgplotbands(S_tilde_ts_Heston, 
                        main="Risk-Neutral Paths - Heston")

# Sample plots
par(mfrow=c(1,3))
matplot(S_tilde_combined_GBM[, sample(ncol(S_tilde_combined_GBM), 200)], 
        type='l', col=rgb(0,0,1,0.1),
        main="Risk-Neutral Paths (GBM)",
        xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_GBM), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)

matplot(S_tilde_combined_SVJD[, sample(ncol(S_tilde_combined_SVJD), 200)], 
        type='l', col=rgb(0,0,1,0.1),
        main="Risk-Neutral Paths (SVJD)",
        xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_SVJD), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)

matplot(S_tilde_combined_Heston[, sample(ncol(S_tilde_combined_Heston), 200)], 
        type='l', col=rgb(0,0,1,0.1),
        main="Risk-Neutral Paths (Heston)",
        xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_Heston), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)

par(mfrow=c(1,1))
# =============================================================================
# 6 - VERIFY RISK-NEUTRAL PROPERTY
# =============================================================================

cat("\n=== RISK-NEUTRAL VERIFICATION ===\n\n")
## 
## === RISK-NEUTRAL VERIFICATION ===
terminal_prices_rn_GBM <- S_tilde_combined_GBM[n_periods, ]
terminal_prices_rn_SVJD <- S_tilde_combined_SVJD[n_periods, ]
terminal_prices_rn_Heston <- S_tilde_combined_Heston[n_periods, ]
capitalized_stock_price <- S0 * exp(r * maturity)

cat("GBM Risk-Neutral Verification:\n")
## GBM Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_GBM), "\n")
## Empirical mean: 127.8758
cat("Difference:", mean(terminal_prices_rn_GBM) - capitalized_stock_price, "\n")
## Difference: -0.5267259
print(t.test(terminal_prices_rn_GBM - capitalized_stock_price))
## 
##  One Sample t-test
## 
## data:  terminal_prices_rn_GBM - capitalized_stock_price
## t = -1.825, df = 1999, p-value = 0.06816
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.09275998  0.03930816
## sample estimates:
##  mean of x 
## -0.5267259
cat("\nSVJD Risk-Neutral Verification:\n")
## 
## SVJD Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_SVJD), "\n")
## Empirical mean: 128.2571
cat("Difference:", mean(terminal_prices_rn_SVJD) - capitalized_stock_price, "\n")
## Difference: -0.1454399
print(t.test(terminal_prices_rn_SVJD - capitalized_stock_price))
## 
##  One Sample t-test
## 
## data:  terminal_prices_rn_SVJD - capitalized_stock_price
## t = -0.35471, df = 1999, p-value = 0.7228
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.9495529  0.6586730
## sample estimates:
##  mean of x 
## -0.1454399
cat("\nHeston Risk-Neutral Verification:\n")
## 
## Heston Risk-Neutral Verification:
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
## Expected terminal price (Q): 128.4025
cat("Empirical mean:", mean(terminal_prices_rn_Heston), "\n")
## Empirical mean: 127.9524
cat("Difference:", mean(terminal_prices_rn_Heston) - capitalized_stock_price, "\n")
## Difference: -0.4501411
print(t.test(terminal_prices_rn_Heston - capitalized_stock_price))
## 
##  One Sample t-test
## 
## data:  terminal_prices_rn_Heston - capitalized_stock_price
## t = -1.3619, df = 1999, p-value = 0.1734
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.0983259  0.1980436
## sample estimates:
##  mean of x 
## -0.4501411
# Visualization comparison
par(mfrow=c(3, 2))
hist(terminal_prices_physical_GBM, breaks=30, col=rgb(1,0,0,0.5),
     main="Terminal Prices: Physical (GBM)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_GBM), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)

hist(terminal_prices_rn_GBM, breaks=30, col=rgb(0,0,1,0.5),
     main="Terminal Prices: Risk-Neutral (GBM)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_GBM), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)

hist(terminal_prices_physical_SVJD, breaks=30, col=rgb(1,0,0,0.5),
     main="Terminal Prices: Physical (SVJD)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_SVJD), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)

hist(terminal_prices_rn_SVJD, breaks=30, col=rgb(0,0,1,0.5),
     main="Terminal Prices: Risk-Neutral (SVJD)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_SVJD), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)

hist(terminal_prices_physical_Heston, breaks=30, col=rgb(1,0,0,0.5),
     main="Terminal Prices: Physical (Heston)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_Heston), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)

hist(terminal_prices_rn_Heston, breaks=30, col=rgb(0,0,1,0.5),
     main="Terminal Prices: Risk-Neutral (Heston)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_Heston), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)

par(mfrow=c(1,1))
# =============================================================================
# 7 - OPTION PRICING
# =============================================================================

cat("\n=== OPTION PRICING ===\n\n")
## 
## === OPTION PRICING ===
bs_price <- function(S, K, r, sigma, T, q = 0) {

  d1 <- (log(S / K) + (r - q + 0.5 * sigma^2) * T) / (sigma * sqrt(T))
  d2 <- d1 - sigma * sqrt(T)

  call <- S * exp(-q * T) * pnorm(d1) - K * exp(-r * T) * pnorm(d2)
  put  <- K * exp(-r * T) * pnorm(-d2) - S * exp(-q * T) * pnorm(-d1)

  list(call = call, put = put)
}


strikes <- seq(80, 160, by=10)
d_f <- exp(-r * maturity)

# Function to price options
price_options <- function(terminal_prices, strikes, discount_factor) {
  n_strikes <- length(strikes)
  call_prices <- numeric(n_strikes)
  bs_call_prices <- numeric(n_strikes)
  put_prices <- numeric(n_strikes)
  bs_put_prices <- numeric(n_strikes)
  call_se <- numeric(n_strikes)
  put_se <- numeric(n_strikes)
  
  for (k in 1:n_strikes) {
    K <- strikes[k]
    call_payoffs <- pmax(terminal_prices - K, 0)
    call_prices[k] <- mean(call_payoffs) * discount_factor
    call_se[k] <- sd(call_payoffs) / sqrt(length(call_payoffs)) * discount_factor
    
    put_payoffs <- pmax(K - terminal_prices, 0)
    put_prices[k] <- mean(put_payoffs) * discount_factor
    put_se[k] <- sd(put_payoffs) / sqrt(length(put_payoffs)) * discount_factor
    
    bs_prices <- bs_price(S0, K, r, sigma=sigma, T=5, q = 0)
    bs_call_prices[k] <- bs_prices$call
    bs_put_prices[k] <- bs_prices$put
  }
  
  list(call_prices = call_prices, put_prices = put_prices,
       bs_call_prices = bs_call_prices, bs_put_prices = bs_put_prices,
       call_se = call_se, put_se = put_se)
}

# Price options for all models
options_GBM <- price_options(terminal_prices_rn_GBM, strikes, d_f)
options_SVJD <- price_options(terminal_prices_rn_SVJD, strikes, d_f)
options_Heston <- price_options(terminal_prices_rn_Heston, strikes, d_f)
kableExtra::kable(as.data.frame(options_GBM))
call_prices put_prices bs_call_prices bs_put_prices call_se put_se
37.2870565 0.0013337 37.6959374 0.0000001 0.2246654 0.0013337
29.5063372 0.0086222 29.9079897 0.0000602 0.2241400 0.0057295
21.7615753 0.0518681 22.1260239 0.0061022 0.2217762 0.0145639
14.2625991 0.3408998 14.4726589 0.1407450 0.2104212 0.0372080
7.7729434 1.6392520 7.6644048 1.1204988 0.1759057 0.0826745
3.2759164 4.9302327 3.0014135 4.2455153 0.1193116 0.1418863
0.9644973 10.4068215 0.8280730 9.8601826 0.0623878 0.1912890
0.1787918 17.4091238 0.1608817 16.9809992 0.0242475 0.2163885
0.0157002 25.0340400 0.0226357 24.6307609 0.0059463 0.2238247
kableExtra::kable(as.data.frame(options_Heston))
call_prices put_prices bs_call_prices bs_put_prices call_se put_se
37.3495392 0.0041721 37.6959374 0.0000001 0.2570826 0.0030129
29.5933710 0.0360118 29.9079897 0.0000602 0.2550569 0.0116692
21.9374713 0.1681198 22.1260239 0.0061022 0.2485708 0.0279110
14.6507009 0.6693573 14.4726589 0.1407450 0.2302646 0.0585080
8.3516001 2.1582644 7.6644048 1.1204988 0.1917575 0.1070118
3.7789938 5.3736659 3.0014135 4.2455153 0.1363246 0.1653930
1.2901727 10.6728526 0.8280730 9.8601826 0.0804251 0.2145030
0.3303330 17.5010208 0.1608817 16.9809992 0.0398569 0.2426599
0.0673177 25.0260132 0.0226357 24.6307609 0.0186772 0.2534207
kableExtra::kable(as.data.frame(options_SVJD))
call_prices put_prices bs_call_prices bs_put_prices call_se put_se
37.7713083 0.1886397 37.6959374 0.0000001 0.3023061 0.0587409
30.1244079 0.3297471 29.9079897 0.0000602 0.2943756 0.0732991
22.6284991 0.6218461 22.1260239 0.0061022 0.2812424 0.0937678
15.4824587 1.2638136 14.4726589 0.1407450 0.2592663 0.1231765
9.2623385 2.8317012 7.6644048 1.1204988 0.2205829 0.1645298
4.6409015 5.9982720 3.0014135 4.2455153 0.1668160 0.2151494
1.8977393 11.0431176 0.8280730 9.8601826 0.1118326 0.2617140
0.6901214 17.6235076 0.1608817 16.9809992 0.0688966 0.2916385
0.2333527 24.9547468 0.0226357 24.6307609 0.0402776 0.3074406